#-------------------------------------------------------------------------------
"This function simulates the model and computes moments to compute Tables 5-9"
#-------------------------------------------------------------------------------
#= Output of the function:
target_mms    = [E[b/y]; E[def]; E[SP]; sd(SP); B/(B+b); Elast]                                           #TABLE 5
untgt_mms     = [sd(C)/sd(Y); sd(TB/Y)/sd(Y); cor(C,Y); cor(TB/Y,Y); cor(SP,Y)]                           #TABLE 6
piBE_mms      = [E[π]; sd(π); sd(dBE)_EVENTS; cor(dπ,dlnSP); cor(dπ,dBE); cor(π,Y)                        #TABLE 7
RP_mms        = [E(RP);  E(RP/SP); R(RP/SP, low y); cor(RP,Y); cor(RP/SP,Y);  std(RP); std(SP, high Z)]   #TABLE 8
Moments_Ctype = [E(b/y, Ctype);  E(SP, Ctype); std(SP, Ctype); E(RP, Ctype); std(RP, Ctype)]              #TABLE 9
#------------------------------------------------------------------------------- =#

#-------------------------------------------------------------------------------
function Simulations_Wrapper(ce::CE_Economy_M; T=300_000, Initial_Type=1,
                             do_plots=true, WINDOWS=true, Compute_η=true)
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# RUN THE SIMULATION
#-------------------------------------------------------------------------------
y_sim, type_sim, b_sim, π_sim, def_sim, def_pol, def_id, q_sim, SP_sim, SP_HR_sim, c_sim, tb_sim, m_sim,m_sim_ctfl, ζ_sim, dlnSP_sim, dBE_sim, b_sim_tmp, inf_bp,
elasticity, ΔSP_sim,ΔBE_sim,Δπ_pol = simul_econ(ce,T,1, T0=5_000, Initial_Type=Initial_Type);
RP_sim = SP_sim - SP_HR_sim;

#-------------------------------------------------------------------------------
#Compute Unconditional Moments
#-------------------------------------------------------------------------------
T0 = 50_000
T1 = 250_000
#Create variables to fill in
def_fre,SP_mean, SP_std,  by_mean, Bv_mean,                                        #Table 5 - Targeted moments
_std_c_gdp,_std_tb_gdp,_cor_c_gdp, _cor_tb_gdp,_cor_SP_gdp,                        #Table 6
_mean_π,_std_π, _cor_π_gdp,                                                        #Table 7 [panel A]
_mean_RP,_mean_RPSP,_std_RP,_std_SP_zH,_mean_RPSP_Low_y, _cor_RP_gdp,_cor_RPSP_gdp =  #Table 8 
0.,0.,0.,0.,0.,   0.,0.,0.,0.,0.,    0.,0.,0.,    0.,0.,0.,0.,0.,0.,0.;


if WINDOWS==true #Targeted & Untargeted Moments (around no-default periods)
def_fre,SP_mean, SP_std,  by_mean, Bv_mean,                                        #Table 5 - Targeted moments
_std_c_gdp,_std_tb_gdp,_cor_c_gdp, _cor_tb_gdp,_cor_SP_gdp,                        #Table 6
_mean_π,_std_π, _cor_π_gdp,                                                        #Table 7 [panel A]
_mean_RP,_mean_RPSP,_std_RP,_std_SP_zH,_mean_RPSP_Low_y, _cor_RP_gdp,_cor_RPSP_gdp =  #Table 8  
compute_moments(ce, y_sim[T0:T1], type_sim[T0:T1], b_sim[T0:T1], π_sim[T0:T1], def_sim[T0:T1],
                    def_pol[T0:T1], def_id[T0:T1], q_sim[T0:T1], SP_sim[T0:T1], SP_HR_sim[T0:T1],
                    c_sim[T0:T1], tb_sim[T0:T1], m_sim[T0:T1],   ζ_sim[T0:T1],
                    dlnSP_sim[T0:T1], dBE_sim[T0:T1], b_sim_tmp[T0:T1], inf_bp[T0:T1],
                    do_plots=do_plots, WINDOWS=WINDOWS);

elseif WINDOWS==false #Targeted moments only (Table 5)
def_fre,SP_mean, SP_std,  by_mean, Bv_mean =
compute_moments(ce, y_sim[T0:T1], type_sim[T0:T1], b_sim[T0:T1], π_sim[T0:T1], def_sim[T0:T1],
                    def_pol[T0:T1], def_id[T0:T1], q_sim[T0:T1], SP_sim[T0:T1], SP_HR_sim[T0:T1],
                    c_sim[T0:T1], tb_sim[T0:T1], m_sim[T0:T1],   ζ_sim[T0:T1],
                    dlnSP_sim[T0:T1], dBE_sim[T0:T1], b_sim_tmp[T0:T1], inf_bp[T0:T1],
                    do_plots=do_plots, WINDOWS=WINDOWS);
end

# -------------------------------------------------------------------------------
# Spread decomposition for the C-type (Table 9)
# -------------------------------------------------------------------------------
Moments_Ctype  = zeros(5,2)
std_y          = std(y_sim[(def_sim.==0)][T0:T1])
mean_y         = mean(y_sim[(def_sim.==0)][T0:T1])
y_low          = mean_y-std_y   
filter_obs     = []
for ii=1:2
  if ii==1
  filter_obs      = (def_sim.==0).*(type_sim.==1)                 #All periods for C-type
  elseif ii==2
  filter_obs      = (def_sim.==0).*(type_sim.==1).*(y_sim.<y_low) #C-type & Output<Avg-STD
  end
by_mean_C             = mean(b_sim[filter_obs]./y_sim[filter_obs])
SP_mean_C             = mean(SP_sim[filter_obs])
SP_std_C              = std(SP_sim[filter_obs]);
RP_mean_C             = mean(RP_sim[filter_obs])
RP_std_C              = std(RP_sim[filter_obs]);
Moments_Ctype[:,ii]   = [ by_mean_C;  SP_mean_C; SP_std_C; RP_mean_C; RP_std_C]*100 ;
end


#-------------------------------------------------------------------------------
# Compute elasticity η(SP,BE)
#-------------------------------------------------------------------------------
elastAVG, sd_dBE_event, cor_dπ_dlnSP, cor_dπ_dBE = 0.,0.,0.,0.;
if Compute_η == true
#Compute moments needed for the GIRF
π_std      = std(π_sim[ (type_sim.==2) .* (def_sim.==0)])
ζ_avg_C    = mean(ζ_sim[(type_sim.==1) .* (def_sim.==0)])
ζ_std_C    = std(ζ_sim[ (type_sim.==1) .* (def_sim.==0)])
ζC_l       = ζ_avg_C-ζ_std_C
ζC_h       = ζ_avg_C+ζ_std_C
#Run the simulation
y_sim, type_sim, b_sim, π_sim, def_sim, def_pol, def_id, q_sim, SP_sim, SP_HR_sim, c_sim, tb_sim, m_sim, m_sim_ctfl, ζ_sim, dlnSP_sim, dBE_sim, b_sim_tmp, inf_bp,
elasticity, ΔSP_sim,ΔBE_sim,Δπ_pol = simul_econ(ce,T,1, T0=5_000, π_shock=π_std)
#Filter observations
filter_sp     = filter_fx(def_sim,type_sim,ζ_sim, SP_sim, ζ_min=ζC_l, ζ_max=ζC_h)
#Compute average elasticity and other moments
elastAVG      = mean(elasticity[filter_sp[:]])
sd_dBE_event  = std(ΔBE_sim[:,1,2][filter_sp[:]])
cor_dπ_dlnSP  = cor(Δπ_pol[filter_sp[:]],ΔSP_sim[:,1,2][filter_sp[:]])
cor_dπ_dBE    = cor(Δπ_pol[filter_sp[:]],ΔBE_sim[:,1,2][filter_sp[:]])
end
#-------------------------------------------------------------------------------


#-------------------------------------------------------------------------------
# Store results
#-------------------------------------------------------------------------------
#Table 5
target_mms   = [by_mean;      def_fre;      SP_mean;         SP_std;           Bv_mean*100;        elastAVG]              # E[b/y]; E[def]; E[SP]; sd(SP); B/(B+b); Elast  [excludes only default periods]
#Table 6
untgt_mms    = [_std_c_gdp;   _std_tb_gdp;  _cor_c_gdp*100;  _cor_tb_gdp*100;  _cor_SP_gdp*100]                           # sd(C)/sd(Y); sd(TB/Y)/sd(Y); c(C,Y); c(TB/Y,Y); c(SP,Y)
#Table 7 [panels A and B]
piBE_mms     = [_mean_π;      _std_π;       sd_dBE_event;    cor_dπ_dlnSP*100; cor_dπ_dBE*100;     _cor_π_gdp*100]        # E[π]; sd(π); sd(dBE)_EVENTS; c(dπ,dlnSP); c(dπ,dBE); c(π,Y)
#Table 8
RP_mms       = [_mean_RP*100;  _mean_RPSP; _mean_RPSP_Low_y; _cor_RP_gdp*100;  _cor_RPSP_gdp*100;  _std_RP*100; _std_SP_zH*100]

return target_mms, untgt_mms, piBE_mms, RP_mms, Moments_Ctype[:,1]
end
#-------------------------------------------------------------------------------
